Only submissions in html format are accepted. Replace NetID by your NetID in the cell below and run the cell. DO NOT CHANGE THE NAME OF THIS NOTEBOOK. If characters have been added during download, rename the file to its original name. Do not modify the last cell, and don't forget to run it. Thanks!

student_id = "solution"
print(student_id)
assignment_name = 'HW4_'+student_id

solution

The following cell should always be the first coding cell of your python notebooks

"""
importing the necessary libraries, do not modify
"""
%matplotlib inline 
# plots graphs within the notebook
%config InlineBackend.figure_format='svg' # not sure what this does, may be default images to svg format

from IPython.display import display,Image, Latex, Math
from __future__ import division
from sympy.interactive import printing
printing.init_printing(use_latex='mathjax')

import time
from IPython.display import clear_output

import SchemDraw as schem
import SchemDraw.elements as e

import matplotlib.pyplot as plt
import numpy as np
import math
import scipy.constants as sc

import sympy as sym

from IPython.core.display import HTML
def header(text):
    raw_html = '<h4>' + str(text) + '</h4>'
    return raw_html

def box(text):
    raw_html = '<div style="border:1px dotted black;padding:2em;">'+str(text)+'</div>'
    return HTML(raw_html)

def nobox(text):
    raw_html = '<p>'+str(text)+'</p>'
    return HTML(raw_html)

def addContent(raw_html):
    global htmlContent
    htmlContent += raw_html
    
class PDF(object):
  def __init__(self, pdf, size=(200,200)):
    self.pdf = pdf
    self.size = size

  def _repr_html_(self):
    return '<iframe src={0} width={1[0]} height={1[1]}></iframe>'.format(self.pdf, self.size)

  def _repr_latex_(self):
    return r'\includegraphics[width=1.0\textwidth]{{{0}}}'.format(self.pdf)

class ListTable(list):
    """ Overridden list class which takes a 2-dimensional list of 
        the form [[1,2,3],[4,5,6]], and renders an HTML Table in 
        IPython Notebook. """
    
    def _repr_html_(self):
        html = ["<table>"]
        for row in self:
            html.append("<tr>")
            
            for col in row:
                html.append("<td>{0}</td>".format(col))
            
            html.append("</tr>")
        html.append("</table>")
        return ''.join(html)
    
font = {'family' : 'serif',
        'color'  : 'black',
        'weight' : 'normal',
        'size'   : 18,
        }

from scipy.constants.constants import C2K
from scipy.constants.constants import K2C
from scipy.constants.constants import F2K
from scipy.constants.constants import K2F
from scipy.constants.constants import C2F
from scipy.constants.constants import F2C

time_stamp = np.zeros(10)
time_stamp[0] = time.time()

Error: Unable to load license data from: '/Users/dubief/.continuum/license_mkl_1458597647.txt'

time_stamp[1] = time.time()

Parameters

The following code solves the 1D pure-conduction heat equation: $$ \frac{\partial T}{\partial t}=\alpha\frac{\partial^2 T}{\partial x^2} $$

over $0\leq x\leq L$ with a time-dependent temperature condition on one side and an isothermal boundary condition on the other: $$ T(x=0,t)=T_m + A\sin(\omega t),\;T(x=L,t)=T_L $$

where, initially, $T_m=300\text{K}=T_L$ and $A=20\text{K}$. The period of oscillation is $t_p = 60\text{s}$, and $\omega=2\pi/t_p$. The length is $2\text{cm}$ and the material is wood with $\rho=545\text{kg/m$^3$}$, $Cp = 2385. \text{J/kg.K}$ and $k = 0.17 \text{W/m.K}$.

The time advancement is explicit Euler and the second spatial derivative is approximated with a second order central scheme: $$ \frac{T_i^{n+1}-T_i^n}{\Delta t}=\alpha\frac{T_{i-1}^n-2T_i^n+T_{i+1}^n}{\Delta x^2} $$

where $\Delta t$ and $\Delta x$ are the time step and the grid step, respectively. The stability of the time advancement scheme requires: $$ \frac{\alpha\Delta t}{\Delta x^2}\leq \frac{1}{2} $$

from Libraries import HT_conduction_extended_surfaces as extsurf
from Libraries import HT_thermal_resistance as res

L = 0.02 #m
T_L = 300. #K
T_m = 300. #K
A = 20. #K
Period = 60 #s
omega = 2*np.pi/Period
rho = 545. #kg/m^3
Cp = 2385. #J/kg.K
k = 0.17 #W/m.K

alpha = k/(rho*Cp)

# number of nodes between x = 0 and x = L

N = 100
x = np.linspace(0.,L,N) #uniform distribution of points
dx = x[1] - x[0]

# time step
dt = 0.01 #s
dt = np.minimum(dt,0.5*dx**2/alpha)

#simulation time
t_max = 2*Period

Ndt = np.int(t_max/dt)


#Allocation of memory

T = np.zeros(N)
T_old = np.zeros(N)

# Initial conditions
T[:] = T_m

# simulation 

t = 0. #time variable
t_plot = 1. # plot temperature distribution every t_plot seconds
Ntplot = np.int(t_max/t_plot)+1
Tplot = np.zeros((Ntplot,N))
tplot = np.zeros(Ntplot)
t_p = 0.
iplot = 0
Tplot[iplot,:] = T[:]
while (t < t_max):
    T_old[:] = np.copy(T)
    t += dt
    t_p += dt
    T[0] = T_m + A*np.sin(omega*t)
    for i in range(1,N-1):
        T[i] = T_old[i] + alpha*dt/dx**2*\
                    (T_old[i-1]-2.*T_old[i]+T_old[i+1])
    T[N-1] = T_L
    
    if (t_p >= t_plot):
        tplot[iplot+1] = t
        Tplot[iplot+1,:] = T[:]
        iplot +=1
        #print(t)
        #plt.title(r"t=%4.0f s" %t, fontdict = font)
        #plt.plot(x,K2F(T),lw=2)
        #plt.ylim(K2F(T_L-A),K2F(T_L+A))
        #plt.xlim(0.,L)
        #plt.grid(True)
        #plt.xlabel(r'$x$', fontdict = font)
        #plt.ylabel(r'$T$ (F)', fontdict = font)
        #plt.show()
        #clear_output(wait=True)
        t_p = 0.

plt.title(r"t=%4.0f s" %t, fontdict = font)
plt.plot(x,K2F(T),lw=2)
plt.ylim(K2F(T_L-A),K2F(T_L+A))
plt.xlim(0.,L)
plt.grid(True)
plt.xlabel(r'$x$', fontdict = font)
plt.ylabel(r'$T$ (F)', fontdict = font)
plt.show()

# JSAnimation import available at https://github.com/jakevdp/JSAnimation
from JSAnimation import IPython_display
from matplotlib import animation

fig = plt.figure(figsize=(6,6))
ax = plt.axes(xlim=(0, L), ylim=(T_L-A, T_L+A))
time_template = 'time = %3.0fs'
time_text = ax.text(0.015, 285, '')
ax.set_xlabel(r'$x$',fontdict= font)
ax.set_ylabel(r'$T$ (K)',fontdict= font)
line, = ax.plot([], [], lw=2)


def init():
    line.set_data([], [])
    time_text.set_text('')
    return line,time_text

def animate(i):
    line.set_data(x, Tplot[i,:])
    time_text.set_text(time_template % (tplot[i]))
    return line,time_text

animation.FuncAnimation(fig, animate, init_func=init,
                        frames=Ntplot, interval=50, blit=True)

Question 1

Run the simulation for $T_L=280\text{K}$. At steady state, i.e. if $A=0$, how should the temperature vary? Could you devise a plot to assess if there is any region of the wall where the temperature at the end of the two periods reaches steady-state like behavior? (Hint: the numpy function np.gradient used to create the variable delta_x below could be useful http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.gradient.html)

from Libraries import HT_conduction_extended_surfaces as extsurf
from Libraries import HT_thermal_resistance as res

L = 0.02 #m
T_L = 280. #K
T_m = 300. #K
A = 20. #K
Period = 60 #s
omega = 2*np.pi/Period
rho = 545. #kg/m^3
Cp = 2385. #J/kg.K
k = 0.17 #W/m.K

alpha = k/(rho*Cp)

# number of nodes between x = 0 and x = L

N = 100
x = np.linspace(0.,L,N) #uniform distribution of points
dx = x[1] - x[0]
delta_x = np.gradient(x)

# time step
dt = 0.001 #s
dt = np.minimum(dt,0.5*dx**2/alpha)

#simulation time
t_max = 20*Period




#Allocation of memory

T = np.zeros(N)
T_steady = np.zeros(N)
T_steady = T_m + (T_L-T_m)*x/L
T_old = np.zeros(N)

# Initial conditions
T[:] = T_m

# simulation 

t = 0. #time variable
t_plot = 10. # plot temperature distribution every t_plot seconds
Ntplot = np.int(t_max/t_plot)+1
Tplot = np.zeros((Ntplot,N))
tplot = np.zeros(Ntplot)
t_p = 0.
iplot = 0
Tplot[iplot,:] = T[:]
def RHS(T):
    return alpha/dx**2*(T_old[0:N-2]-2.*T_old[1:N-1]+T_old[2:N])
while (t < t_max):
    T_old[:] = np.copy(T)
    t += dt
    t_p += dt
    
    T[1:N-1] = T_old[1:N-1] + dt*RHS(T_old)
    
    T[0] = T_m + A*np.sin(omega*t)
    T[N-1] = T_L
    
    if (t_p >= t_plot):
        #print(t)
        tplot[iplot+1] = t
        Tplot[iplot+1,:] = T[:]
        iplot +=1
        #plt.title(r"t=%4.0f s" %t, fontdict = font)
        #plt.plot(x,K2F(T),lw=2)
        #plt.plot(x,K2F(T_steady),lw=2)
        #plt.ylim(K2F(T_m-A),K2F(T_m+A))
        #plt.xlim(0.,L)
        #plt.grid(True)
        #plt.xlabel(r'$x$', fontdict = font)
        #plt.ylabel(r'$T$ (F)', fontdict = font)
        #plt.show()
        #clear_output(wait=True)
        t_p = 0.


plt.title(r"t=%4.0f s" %t, fontdict = font)
plt.plot(x,K2F(T),lw=2)
plt.plot(x,K2F(T_steady),lw=2)
plt.ylim(K2F(T_m-A),K2F(T_m+A))
plt.xlim(0.,L)
plt.grid(True)
plt.xlabel(r'$x$', fontdict = font)
plt.ylabel(r'$T$ (F)', fontdict = font)
plt.show()

# JSAnimation import available at https://github.com/jakevdp/JSAnimation
from JSAnimation import IPython_display
from matplotlib import animation

fig = plt.figure(figsize=(6,6))
ax = plt.axes(xlim=(0, L), ylim=(T_m-A, T_m+A))
time_template = 'time = %4.0fs'
time_text = ax.text(0.015, 315, '')
ax.set_xlabel(r'$x$',fontdict= font)
ax.set_ylabel(r'$T$ (K)',fontdict= font)
steady = ax.plot(x,T_steady, 'r-', lw=2)
line, = ax.plot([], [], 'b-', lw=2)


def init():
    line.set_data([], [])
    time_text.set_text('')
    return line,time_text

def animate(i):
    line.set_data(x, Tplot[i,:])
    time_text.set_text(time_template % (tplot[i]))
    return line,time_text

animation.FuncAnimation(fig, animate, init_func=init,
                        frames=Ntplot, interval=50, blit=True)

""" Plot to be used to assess the existence of steady-state like behavior 
in the temperature distribution at the end of two periods"""
plt.plot(x,T,lw=2)
plt.plot(x,(T_steady),lw=2)
plt.xlim(0.0,L)
plt.grid(True)
plt.xlabel(r'$x$', fontdict = font)
plt.ylabel(r'$T$ (K)', fontdict = font)
plt.show()

Question 2

Compute the heat flux at $x=0$ and $x=L$.

#Here calculate and print out the heat fluxes at x=0 and x=L
q_0 = -k * (T[1] - T[0])/dx
print("Heat flux at x=0 %3.2f W" %q_0)
q_L = -k * (T[N-1] - T[N-2])/dx
print("Heat flux at x=L %3.2f W" %q_L)

Heat flux at x=0 2048.08 W
Heat flux at x=L 177.31 W

Question 3

Change the boundary condition at $x=L$ to a convection boundary condition with $h=100\text{W/m$^2$.K}$ and $T_\infty=280\text{K}$.

Question 4

Run the simulation for $N=10,20,50,100,500,1000$ and calculate the relative error in both heat fluxes (at each extremity) of the smaller grids relative to the the $N=1000$-grid. (What to do: for each simulations, compute the fluxes and fill the arrays q_0_tab and q_L_tab in the cell below.

from Libraries import HT_conduction_extended_surfaces as extsurf
from Libraries import HT_thermal_resistance as res

L = 0.02 #m
T_m = 300. #K
T_infty = 280. #K
A = 20. #K
Period = 60 #s
omega = 2*np.pi/Period
rho = 545. #kg/m^3
Cp = 2385. #J/kg.K
k = 0.17 #W/m.K
h = 100. #W/m^2.K


alpha = k/(rho*Cp)

# number of nodes between x = 0 and x = L

N = 100
x = np.linspace(0.,L,N) #uniform distribution of points
dx = x[1] - x[0]

# time step
dt = 0.01 #s
dt = np.minimum(dt,0.5*dx**2/alpha)

#simulation time
t_max = 2*Period




#Allocation of memory

T = np.zeros(N)
T_old = np.zeros(N)

# Initial conditions
T[:] = T_m

# simulation 

t = 0. #time variable
t_plot = 1. # plot temperature distribution every t_plot seconds
Ntplot = np.int(t_max/t_plot)+1
Tplot = np.zeros((Ntplot,N))
tplot = np.zeros(Ntplot)
t_p = 0.
iplot = 0
Tplot[iplot,:] = T[:]

def RHS(T):
    return alpha/dx**2*(T[0:N-2]-2.*T[1:N-1]+T[2:N])
while (t < t_max):
    T_old[:] = np.copy(T)
    t += dt
    t_p += dt
    
    T[1:N-1] = T_old[1:N-1] + dt * RHS(T_old)
    # Boundary conditions
    T[0] = T_m + A*np.sin(omega*t)
    T[N-1] = (h*T_infty+k*T[N-2]/dx)/(h+k/dx)
    
    if (t_p >= t_plot):
        #print(t)
        tplot[iplot+1] = t
        Tplot[iplot+1,:] = T[:]
        iplot +=1
        #plt.title(r"t=%4.0f s" %t, fontdict = font)
        #plt.plot(x,K2F(T),lw=2)
        #plt.plot(x,K2F(T_steady),lw=2)
        #plt.ylim(K2F(T_m-A),K2F(T_m+A))
        #plt.xlim(0.,L)
        #plt.grid(True)
        #plt.xlabel(r'$x$', fontdict = font)
        #plt.ylabel(r'$T$ (F)', fontdict = font)
        #plt.show()
        #clear_output(wait=True)
        t_p = 0.


plt.title(r"t=%4.0f s" %t, fontdict = font)
plt.plot(x,K2F(T),lw=2)
plt.ylim(K2F(T_m-A),K2F(T_m+A))
plt.xlim(0.,L)
plt.grid(True)
plt.xlabel(r'$x$', fontdict = font)
plt.ylabel(r'$T$ (F)', fontdict = font)
plt.show()

# JSAnimation import available at https://github.com/jakevdp/JSAnimation
from JSAnimation import IPython_display
from matplotlib import animation

fig = plt.figure(figsize=(6,6))
ax = plt.axes(xlim=(0, L), ylim=(T_m-A, T_m+A))
time_template = 'time = %3.0fs'
time_text = ax.text(0.015, 285, '')
ax.set_xlabel(r'$x$',fontdict= font)
ax.set_ylabel(r'$T$ (K)',fontdict= font)
line, = ax.plot([], [], lw=2)


def init():
    line.set_data([], [])
    time_text.set_text('')
    return line,time_text

def animate(i):
    line.set_data(x, Tplot[i,:])
    time_text.set_text(time_template % (tplot[i]))
    return line,time_text

animation.FuncAnimation(fig, animate, init_func=init,
                        frames=Ntplot, interval=50, blit=True)

#Here calculate and print out the heat fluxes at x=0 and x=L
q_0 = -k * (T[1] - T[0])/dx
print("Heat flux at x=0 %3.2f W" %q_0)
q_L = -k * (T[N-1] - T[N-2])/dx
print("Heat flux at x=L %3.2f W" %q_L)

Heat flux at x=0 1867.63 W
Heat flux at x=L 445.25 W

"""Cell where you compute the relative error"""
q_0_N1000 = 2106.10
q_L_N1000 = 448.39
print(dx,L/(500.-1.))
delta_tab = np.array([L/(10.-1.),L/(20.-1.),L/(100.-1.),L/(500.-1.)])
q_0_tab = np.array([288.73,951.22,1867.13,2079.06])
q_L_tab = np.array([410.46,430.65,445.25,448.04])

delta_one = 0.1*np.power(delta_tab,1)/delta_tab[3]**1
delta_square = 0.1*np.power(delta_tab,2)/delta_tab[3]**2

error_0 = np.abs(q_0_tab-q_0_N1000)/q_0_N1000*100.
error_L = np.abs(q_L_tab-q_L_N1000)/q_L_N1000*100.

plt.loglog(delta_tab,error_0,lw=2, label=r"$x=0$")
plt.loglog(delta_tab,error_L,lw=2, label=r"$x=L$")
plt.loglog(delta_tab,delta_one, 'k:',lw=2, label=r"$\propto\Delta_x$")
plt.loglog(delta_tab,delta_square, 'k--',lw=2, label=r"$\propto\Delta_x^2$")
#plt.ylim(xx,yy)
#plt.xlim(0.0,L)
plt.grid(True)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=2, shadow=True, fancybox=True)
plt.xlabel(r'$x$', fontdict = font)
plt.ylabel(r'relative error (%)', fontdict = font)
plt.show()

0.00020202020202 4.0080160320641287e-05

Question 5

Modify the time advancement scheme to the second order Runge-Kutta, with an error of order magnitude $\Delta t^2$ compared to $\Delta t$ for the explicit Euler: \begin{eqnarray} T_i^{n+1/2} &=& T_i^n+\frac{\Delta t}{2}\left(RHS_i^n\right)\ T_i^{n+1} &=& T_i^n+\Delta t\left(RHS_i^{n+1/2}\right) \end{eqnarray}

Choose $N=100$ and compare the solution with the explicit Euler solution for the same $N$.

from Libraries import HT_conduction_extended_surfaces as extsurf
from Libraries import HT_thermal_resistance as res

L = 0.02 #m
T_m = 300. #K
A = 20. #K
Period = 60 #s
omega = 2*np.pi/Period
rho = 545. #kg/m^3
Cp = 2385. #J/kg.K
k = 0.17 #W/m.K
h = 100. #W/m^2.K
T_infty = 280.

alpha = k/(rho*Cp)

# number of nodes between x = 0 and x = L

N = 100
x = np.linspace(0.,L,N) #uniform distribution of points
dx = x[1] - x[0]

# time step
dt = 0.001 #s
dt = np.minimum(dt,0.5*dx**2/alpha)

#simulation time
t_max = 2*Period




#Allocation of memory

T = np.zeros(N)
T_old = np.zeros(N)

# Initial conditions
T[:] = T_m

# simulation 

t = 0. #time variable
t_plot = 1. # plot temperature distribution every t_plot seconds
Ntplot = np.int(t_max/t_plot)+1
Tplot = np.zeros((Ntplot,N))
tplot = np.zeros(Ntplot)
t_p = 0.
iplot = 0
Tplot[iplot,:] = T[:]

def RHS(T):
    return alpha/dx**2*(T[0:N-2]-2.*T[1:N-1]+T[2:N])

while (t < t_max):
    
    t += dt
    t_p += dt
    T_old[:] = np.copy(T)
    
    T[1:N-1] = T_old[1:N-1] + 0.5*dt*RHS(T_old)
    T[1:N-1] = T_old[1:N-1] + dt*RHS(T)
    
    T[0] = T_m + A*np.sin(omega*t)
    T[N-1] = (h*T_infty+k*T[N-2]/dx)/(h+k/dx)
    
    if (t_p >= t_plot):
        #print(t)
        tplot[iplot+1] = t
        Tplot[iplot+1,:] = T[:]
        iplot +=1
        #plt.title(r"t=%4.0f s" %t, fontdict = font)
        #plt.plot(x,K2F(T),lw=2)
        #plt.plot(x,K2F(T_steady),lw=2)
        #plt.ylim(K2F(T_m-A),K2F(T_m+A))
        #plt.xlim(0.,L)
        #plt.grid(True)
        #plt.xlabel(r'$x$', fontdict = font)
        #plt.ylabel(r'$T$ (F)', fontdict = font)
        #plt.show()
        #clear_output(wait=True)
        t_p = 0.


plt.title(r"t=%4.0f s" %t, fontdict = font)
plt.plot(x,K2F(T),lw=2)
plt.ylim(K2F(T_m-A),K2F(T_m+A))
plt.xlim(0.,L)
plt.grid(True)
plt.xlabel(r'$x$', fontdict = font)
plt.ylabel(r'$T$ (F)', fontdict = font)
plt.show()

#Here calculate and print out the heat fluxes at x=0 and x=L
q_0 = -k * (T[1] - T[0])/dx
print("Heat flux at x=0 %3.2f W" %q_0)
q_L = -k * (T[N-1] - T[N-2])/dx
print("Heat flux at x=L %3.2f W" %q_L)

Heat flux at x=0 1867.84 W
Heat flux at x=L 445.26 W

!ipython nbconvert --to html ME144-HW4-solution1.ipynb --template full --output $assignment_name

[NbConvertApp] Converting notebook ME144-TH_1_solution.ipynb to html
[NbConvertApp] Writing 308210 bytes to TH_1_solution.html